#include <TROOT.h>

void Z_CHIRALVAL(char *momentum_order, int beta100, int l, int t, int nclust, char *dir, int imom, char  *ensemble, char *pmprint, char *Oa2print, char *gammaprint , int iprop, int useMPS,int subGP)
{
	
        L=l;
        T=t;
        Nclust=nclust;

	printf("Chiral fit for Z%s ", gammaprint);
	if (useMPS) printf(" using MPS^2");
	else printf(" using mpolar");
	if (subGP) printf(" and substracting 1 GP\n");
	else printf(" \n");

        double ZAPREV=0.757;
        int Nmass;
	int set_p[4];
	char setpprint[16];

         /////////////////////////////
        // INPUT FILES///////////
        ////////////////////////////////



        // mass input files

        char mu_dat[1024];
        sprintf(mu_dat,"../INPUT/mu_%s%s_%d.dat",ensemble,pmprint,beta100);
        FILE *mudat=fopen(mu_dat,"r");
        fscanf(mudat,"%i",&Nmass);
        double *mu=(double*)malloc(sizeof(double)*(Nmass));
        for (int im=0; im<Nmass; im++)
        fscanf(mudat,"%lf",&mu[im]);
		

       char mPCAC_dat[1024];
       sprintf(mPCAC_dat,"../INPUT/mPCAC_%s%s_%d.dat", ensemble,pmprint,beta100);
       double mPCAC;
       FILE *mPCACdat=fopen(mPCAC_dat,"r");
       fscanf(mPCACdat,"%lf",&mPCAC);
       fscanf(mPCACdat,"%lf",&ZAPREV);

       // total non equivalent combination of Nmass: Nmass_comb
       // I checked that Z(m1,m2)=Z(m2,m1) 
       int Nmass_comb=Nmass;
       for (int i=1; i<Nmass; i++)  Nmass_comb=Nmass_comb+(Nmass-i); 

       if (useMPS)
       {
	int aux_int;
	double aux_double;
        char MPS_dat[1024];
        sprintf(MPS_dat,"../INPUT/MPS_%s%s_%d.dat",ensemble,pmprint,beta100);
        FILE *MPSdat=fopen(MPS_dat,"r");

        double **MPS=(double**)malloc(sizeof(double*)*(Nmass));
	for (int im1=0; im1<Nmass; im1++) MPS[im1]=(double*)malloc(sizeof(double)*(Nmass));

	for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++){
        fscanf(MPSdat,"%lf ",&MPS[im1][im2]);
	}
       }
       /////////////////////////
       //// Momenta set
      //////////////////////////////

    whichProp(iprop,setpprint);
    whichSetp(iprop,set_p);
    char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
    read_momenta_data(momentafile);
    double *nn=(double*)malloc(sizeof(double)*(Nmom));
    compute_nn(momentum_order, nn);
    Momentum_Magro(momentum_order, set_p);
    // equivalent momenta
    Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
    Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);



     //////////////////////////////////
     // IO files
     ///////////////////////////////

     char file[1024];
     char dir_qav[1024];

     sprintf(dir_qav,"%s/%s/QAV/",dir,setpprint);
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

  
	//////////////////////////////////
	// Read Zs
	/////////////////////////////
	int Nm;
	if (strcmp(gammaprint,"q")==0) Nm=Nmass;
	else Nm=Nmass_comb;

        double  **Z=malloc(sizeof(double*)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Z[iclust]=malloc(sizeof(double)*(Nm));

	double *aux_clust=malloc(sizeof(double)*(Nclust+1));

	if (strcmp(gammaprint,"q")==0) 	
	{
       		for (int imass=0; imass<Nm; imass++)
    		{
	         sprintf(file,"Z%s_%s_m%d",gammaprint,Oa2print,imass);
        	 read_input_onemomenta(dir_qav,file,aux_clust,imom);
		 for (int iclust=0; iclust<Nclust+1; iclust++) Z[iclust][imass]=aux_clust[iclust];
    		}
	}
	else
	{
		int imass=0;
		for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
        	{
		 sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
        	 read_input_onemomenta(dir_qav,file,aux_clust,imom);
	         for (int iclust=0; iclust<Nclust+1; iclust++) Z[iclust][imass]=aux_clust[iclust];
		imass++;
       		}
    	}


	 //for (int imass=0; imass<Nm; imass++) printf("%8.5f\n",Z[Nclust][imass]);



      // auxiliar for minuit routines
      double *ZAUX=malloc(sizeof(double)*(Nm));
      double **ZAUX_CLUST=malloc(sizeof(double*)*(Nm));
      for (int imass=0; imass<Nm; imass++) ZAUX_CLUST[imass]=malloc(sizeof(double)*(Nclust+1));

      double *MAUX=malloc(sizeof(double)*(Nm));
     // double *MAUXPOL=malloc(sizeof(double)*(Nm));


      if (useMPS)
      {
	  if (strcmp(gammaprint,"q")==0) for (int imass=0; imass<Nmass; imass++) MAUX[imass]=pow(MPS[imass][imass],2);
	  else {
		int imass=0;
                for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++){ 
			MAUX[imass]=pow(MPS[im1][im2],2);
	//		MAUXPOL[imass]=sqrt(pow(mu[im1],2)+pow(ZAPREV*mPCAC,2))+sqrt(pow(mu[im2],2)+pow(ZAPREV*mPCAC,2));
			imass++;
		}
	  }
      }
      if (!useMPS)
      {
      	if (strcmp(gammaprint,"q")==0)  for (int imass=0; imass<Nmass; imass++) MAUX[imass]=sqrt(pow(mu[imass],2)+pow(ZAPREV*mPCAC,2));
	else {
		int imass=0;
		for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++){ 
			MAUX[imass]=sqrt(pow(mu[im1],2)+pow(ZAPREV*mPCAC,2))+sqrt(pow(mu[im2],2)+pow(ZAPREV*mPCAC,2));
			imass++;
			}
	}
      }	
      double  *ZCHI=malloc(sizeof(double)*(Nclust+1)); 
      double  *invZCHI=malloc(sizeof(double)*(Nclust+1));
      double  *ZLIN=malloc(sizeof(double)*(Nclust+1));
      double  *GP=malloc(sizeof(double)*(Nclust+1));
      double  *chi2=malloc(sizeof(double)*(Nclust+1));

	

      sprintf(file,"Z%s_%s",gammaprint,Oa2print);

      for (int iclust=0; iclust<Nclust+1; iclust++)
      {
		if (!subGP) 
		{
		        for (int imass=0; imass<Nm; imass++){
        	                 ZAUX[imass]=Z[iclust][imass];
                	         for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=Z[iclust_loop][imass];
	                }
			lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nm,Nclust,ZCHI[iclust],ZLIN[iclust],chi2[iclust]);
			GP[iclust]=0.;
			write_output_chiral_extrapolation_0GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nm,ZCHI,ZLIN,chi2);
		}
		if (subGP) 
		{
	                for (int imass=0; imass<Nm; imass++){
        	                 ZAUX[imass]=1./Z[iclust][imass];
                	         for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=1./Z[iclust_loop][imass];
                	}
 	                lin_chiralfit_withGP(ZAUX,ZAUX_CLUST,MAUX,Nm,Nclust,invZCHI[iclust],ZLIN[iclust],GP[iclust],chi2[iclust]);
			 // lin_chiralfit_2_withGP(ZAUX,ZAUX_CLUST,MAUX,MAUXPOL, Nm,Nclust,invZCHI[iclust],ZLIN[iclust],GP[iclust],chi2[iclust]);

	                ZCHI[iclust]=1./invZCHI[iclust];
                        write_output_chiral_extrapolation_1GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nm,ZCHI,ZLIN,GP,chi2);
	
		}
    }


   // Output 

   for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZCHI[iclust];
   write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq,imom);
  
    printf("Output for %s written in %s\n",file,dir_qav);


    delete Z;
    delete aux_clust;
    delete ZAUX; delete ZAUX_CLUST; delete MAUX;
    delete ZCHI; delete invZCHI; delete ZLIN; delete GP;
  

}
